function [parms_fitted, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_ver4(parms, kappa_guess, b_guess, eta_guess, z1_guess, chi_guess, a_guess, b_sigma_x_REC, ...
        U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, nom_yld_1y_vol_target, ave_mu3_vol_target, NumTries, fsolve_opts)
% target moments:
% E[U], vol[u], vol(N*w)/vol(Y), vol(dlogY), vol(ylds), vol(mu3_1y)
% parameters:
% kappa, b, eta, z(1), chi, sigma_x_REC
% sigma_x(N,z) = exp(a + sigma_x_REC*1{z==1})

    if z1_guess>=0; error('z1_guess should be <0'); end
    if eta_guess<=0 || eta_guess>=1; error('check eta_guess'); end
    if kappa_guess<=0; error('check kappa0_guess'); end
    if b_guess<=0; error('check b_guess'); end
    if chi_guess>1 || chi_guess<0; error('check chi_guess'); end

    if parms.Nz~=2; error('Code is for 2 states.'); end
    if abs(parms.grid_z(2)-0)>1e-10; error('Code assumes z(2)=0.'); end
    
%     fixedpt_tol = 1e-8;
    % fixedpt_tol = 1e-7;
    fixedpt_tol = 1e-6; % speed things up

    inflation_process.mu_pi = 6.1592e-04;
    inflation_process.rho_pi = 0.81;
    inflation_process.sigma_pi = 0.00245;
    inflation_process.theta_pi = -0.338;
    inflation_process.loading_dlogC = -0.035;

    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;
    [NN_mesh,xx_mesh] = meshgrid(parms.grid_N(:),parms.grid_x(:));
    Pzz_Nxzz = permute(repmat(parms.StateTransitionProbs,[1 1 NN Nx]),[3 4 1 2]);
    mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
    var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi + inflation_process.theta_pi^2)...
        *inflation_process.sigma_pi^2/(1 - inflation_process.rho_pi^2);
    var_e = inflation_process.sigma_pi^2;
    cov_Xe = inflation_process.sigma_pi^2;
    
    N_Nxz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]);
    N_Nz = repmat(parms.grid_N(:),[1 parms.Nz]);
    U_Nz = 1 - N_Nz;
    U_Nxz = 1 - N_Nxz;
    x_Nxz = permute(repmat(parms.grid_x(:),[1 parms.NN parms.Nz]),[2 1 3]);
    
    phi_of_x = @(x) 1./(1 + exp(-x));
    phi_Nxz = phi_of_x(x_Nxz);
    log_phi_Nxz = log(phi_Nxz);
    input_mu3_Nxz = (1 - U_Nxz).*U_Nxz.*(1 - 2*U_Nxz).*(log_phi_Nxz.^3);
    
    FLAG_SUCCESSFUL_FSOLVE = false;
    
    for iter = 1:NumTries
        fprintf('iter %g with kappa_guess=%g, b_guess=%g, eta_guess=%g, z1_guess=%g, chi_guess=%g, a_guess=%g.\n', ...
            iter, kappa_guess, b_guess, eta_guess, z1_guess, chi_guess, a_guess);
        [XFit, FitVal, EXITFLAG] = fsolve(@FitFun,[log(kappa_guess); ...
            log(b_guess); log(1/eta_guess-1); log(-z1_guess); ...
            log(1/chi_guess - 1); a_guess],fsolve_opts);
        if EXITFLAG>0
            fprintf('Found roots on iteration %g.\n', iter);
            FLAG_SUCCESSFUL_FSOLVE = true;
            break;
        else
            fprintf('iter %g, fitted vals (not solution):\n kappa=%g, b=%g, eta=%g, z1=%g, chi=%g, a=%g.\n', ...
                iter, exp(XFit(1)), exp(XFit(2)), 1/(1 + exp(XFit(3))), ...
                -exp(XFit(4)), 1/(1 + exp(XFit(5))), XFit(6));
            kappa_guess = 0.1 + 5*rand;
            b_guess = 0.3 + 0.6*rand;
            eta_guess = 0.1*rand;
            z1_guess = fit_std_dz_2state(parms,vol_gy_target/sqrt(2.5+1*rand));
            chi_guess = 1/(1+4*rand);
            a_guess = log(0.01 + 0.2*rand);
        end
    end
    
    % do a last compute
    parms_fitted = parms;
    parms_fitted.kappa(:) = exp(XFit(1));
    parms_fitted.b(:) = exp(XFit(2));
    parms_fitted.eta(:) = 1/(1 + exp(XFit(3)));
    parms_fitted.grid_z(1) = -exp(XFit(4));
    parms_fitted.chi = 1/(1 + exp(XFit(5)));
    parms_fitted.sigma_x_a = XFit(6);
    parms_fitted.sigma_x_REC = b_sigma_x_REC;
    parms_fitted.sigma_x(:,1) = exp(XFit(6) + b_sigma_x_REC);
    parms_fitted.sigma_x(:,2) = exp(XFit(6));

    fprintf('kappa=%g, b=%g, eta=%g, z1=%g, chi=%g, a=%g.\n', ...
        parms_fitted.kappa(1), parms_fitted.b(1), parms_fitted.eta(1), ...
        parms_fitted.grid_z(1), parms_fitted.chi, XFit(6));

    gesol = SolveFixedPoint_Model_20220926(parms_fitted,fixedpt_tol);
    parms_fitted.FitVal = FitVal;

    if ~FLAG_SUCCESSFUL_FSOLVE
        warning('Model fit: not all moments are satisfied.');
    end

    function Z = FitFun(X)
        
        Z = NaN*zeros(6,1);
        
        parms_temp = parms;
        parms_temp.kappa(:) = exp(X(1));
        parms_temp.b(:) = exp(X(2));
        parms_temp.eta(:) = 1/(1 + exp(X(3)));
        parms_temp.grid_z(1) = -exp(X(4));
        parms_temp.chi = 1/(1 + exp(X(5))); 
        parms_temp.sigma_x(:,1) = exp(X(6) + b_sigma_x_REC);
        parms_temp.sigma_x(:,2) = exp(X(6));
        
        z_Nxz = permute(repmat(parms_temp.grid_z(:),[1 parms.NN parms.Nx]),[2 3 1]); % note: this must be updated when z changes!!
        Y_Nxz = exp(z_Nxz).*N_Nxz;

        gesol_temp = SolveFixedPoint_Model_20220926(parms_temp,fixedpt_tol);
        if gesol_temp.VFIerr>fixedpt_tol
            warning('FitModel20220926_ver1: VFI error = %g is greater than tolerance of %g.', ...
                gesol_temp.VFIerr, fixedpt_tol)
            bComputeMoments = false;
        else
            bComputeMoments = true;
        end
        
        %% Newer version
        
        if bComputeMoments

        [Phi_Nxz,ss_dist_Nxz] = MakeSparseTransMatrix_Nxz(parms_temp,gesol_temp);

        % NEW CODE: quarterly average
        [mean_U, std_U] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,U_Nxz);
        [mean_Y, std_Y] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,Y_Nxz);
        [mean_Nw, std_Nw] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,N_Nxz.*gesol_temp.wages);
        ratio_stdNw_stdY = std_Nw/std_Y;
        
        % quarterly output vol (exact), from library\exact moments
        vol_gy = moment_quarterly_output_growth_vol_Nxz(parms_temp,gesol_temp,Phi_Nxz,ss_dist_Nxz);
        
        % cross-sectional mu3 vol
        std_mu3_cons_1y = CalcTimeSeriesDiffVol(input_mu3_Nxz,ss_dist_Nxz,Phi_Nxz,12);

        % nominal rate
        x_next_Nxzz = gen_law_of_motion_x(parms_temp);
        Nnext_Nxz = gen_law_of_motion_N(parms_temp,gesol_temp);

        C_next = zeros(NN,Nx,Nz,Nz);
        for iznext = 1:parms.Nz
            C_next(:,:,:,iznext) = reshape(interp2(NN_mesh,xx_mesh,gesol_temp.C(:,:,iznext)',...
                Nnext_Nxz(:),x_next_Nxzz(:,iznext)),[NN Nx Nz]);
        end
        inflation_process_temp = inflation_process;
        inflation_process_temp.F_Nxzz = inflation_process_temp.loading_dlogC*(C_next - repmat(sum(Pzz_Nxzz.*C_next,4),[1 1 1 Nz]));
        
        gesol_temp.nom_bonds = PriceNomBonds_Model20220926(parms_temp, gesol_temp, inflation_process_temp, 12);

        NY = 5;
        [Phi_NxYz,grid_Y,prob_NxYz] = MakeSparseTransMatrix_NxYz(parms_temp,gesol_temp,inflation_process_temp.rho_pi,...
            inflation_process_temp.F_Nxzz,NY,ss_dist_Nxz);

        idx = find(gesol_temp.nom_bonds.maturities==12); mat_yr = 1;

        temp_input = permute(repmat(-log(gesol_temp.nom_bonds.P_A(:,:,:,idx))/mat_yr,[1 1 1 NY]),[1 2 4 3]) ...
            + permute(repmat(gesol_temp.nom_bonds.P_C_a(idx)*grid_Y(:)/mat_yr,[1 NN Nx Nz]),[2 3 1 4]);
        [temp_mean,temp_std] = Calc_monthly_vol(prob_NxYz,temp_input);
        var_Xe_component = (gesol_temp.nom_bonds.P_B_b(idx)^2*var_X + gesol_temp.nom_bonds.P_B_c(idx)^2*var_e ...
            + 2*gesol_temp.nom_bonds.P_B_b(idx)*gesol_temp.nom_bonds.P_B_c(idx)*cov_Xe)/(mat_yr^2);

        nom_yld_std = 100*sqrt(temp_std^2 + var_Xe_component);
        
        Z(1) = (mean_U - U_mean_target)/U_mean_target;
        Z(2) = (std_U - U_std_target)/U_std_target;
        Z(3) = (ratio_stdNw_stdY - wagebillratio_target)/wagebillratio_target;
        Z(4) = (vol_gy - vol_gy_target)/vol_gy_target;
        Z(5) = (nom_yld_std - nom_yld_1y_vol_target)/nom_yld_1y_vol_target;
        Z(6) = (std_mu3_cons_1y - ave_mu3_vol_target)/ave_mu3_vol_target;
        
        end
        
    end

end